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This study is aimed at development of a numerical model for a beta-type Stirling engine with rhombic- 
drive mechanism. By taking into account the non-isothermal effects, the effectiveness of the regenerative 
channel, and the thermal resistance of the heating head, the energy equations for the control volumes in 
the expansion chamber, the compression chamber, and the regenerative channel can be derived and 
solved. Meanwhile, a fully developed flow velocity profile in the regenerative channel, in terms of the 
reciprocating velocity of the displacer and the instantaneous pressure difference between the expansion 
and the compression chambers, is derived for calculation of the mass flow rate through the regenerative 
channel. In this manner, the internal irreversibility caused by pressure difference in the two chambers 
and the viscous shear effects due to the motion of the reciprocating displacer on the fluid flow in the 
regenerative channel gap are included. Periodic variation of pressures, volumes, temperatures, masses, 
and heat transfers in the expansion and the compression chambers are predicted. A parametric study of 
the dependence of the power output and thermal efficiency on the geometrical and physical parameters, 
involving regenerative gap, distance between two gears, offset distance from the crank to the center of 
gear, and the heat source temperature, has been performed. 

© 2010 Elsevier Ltd. All rights reserved. 


1. Introduction 

Stirling engines are power machines that operate over a closed, 
regenerative thermodynamic cycle, with cyclic compression and 
expansion of the working fluid at different temperature levels, as 
described by Walker [1 j. The working fluid in the Stirling engines 
might be air, nitrogen, helium, or hydrogen. In theory, ideal Stirling 
engines features high thermal efficiency, low emissions, and 
quietness. Most importantly, a variety of heat sources can be 
utilized, including solar energy, waste heat, and fossil fuels. In 
particular, the solar dish power systems (for example, the SES 25- 
kW SunCatcher [2]) typically consisting a solar dish concentrator, 
a sun tracker, and a Stirling engine have received great attention 
from the energy-related researchers in recent years. 

In general, the Stirling engines are classified into three types of 
configuration [3] as follows: 

(1) a-type — this type of configuration features two pistons, each in 
its own cylinder: 

(2) (3-type — this type of configuration has a piston and a displacer 
in the same cylinder; and 
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(3) y-type — this type of configuration has a piston and a displacer, 
each in its own cylinder. 

In addition, a double-acting Stirling engine has multiple cylin¬ 
ders and elongated power pistons, and can be considered as 
a coupled a engines with thermodynamic cycle taking place 
between the top of one piston and the bottom of the next piston. 

Stirling engine technology has come a long way in the past 
several decades. However, new concepts and designs continue to 
appear, according to the review of Ross [4]. Among all, the 
domestic-scale Stirling electric power generator [5] is particularly 
of great market potential. 

In an ideal Stirling engine, compression and expansion cham¬ 
bers were assumed to be maintained at constant volume or 
isothermal conditions during the heating and the cooling 
processes. However, in 1859 Rankine [6] already mentioned that 
neither the heating nor the cooling takes place exactly at constant 
volumes or at constant temperatures. A simple second-order 
analysis of the ideal Stirling engine cycle was proposed by 
G. Schmidt in 1871 [7]. The Schmidt analysis has been used widely 
as an approximation of Stirling engine performance, and Bercho- 
witz and Urieli [8] gave a detailed description of it. The Schmidt 
analysis takes the isothermal analysis of Stirling engines one step 
further by assuming that the volume of the expansion and the 
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Nomenclature 

Ra 

thermal resistance of cooling jacket (°C/kW) 



t 

time (s) 

Ci,C 2 

coefficient of velocity equation 

fmax 

maximum computation time (s) 

Cv 

constant volume specific heat (kj/kg I<) 

f P 

period (s) 

Cp 

constant pressure specific heat (kj/kg I<) 

£o 

reference time (s) 

dr 

diameter of displacer (m) 

T 

temperature (K) 


core diameter of cylinder (m) 

Th 

heat source temperature (K) 

G 

regenerative channel gap (m) 

T l 

heat sink temperature (K) 

h 

enthalpy (kj/kg) 

U 

internal energy (kj/s) 

/d 

height of displacer (m) 

V 

velocity (m/s) 

k 

height of the engine (m) 

V 

volume (m 3 ) 

i|. h, h: 

f 4 lengths of the linkages of the rhombic-drive 

W 0ut 

power output (kj/s) 


mechanism (m) 

w out 

network output per cycle (kj/cycle) 

L 

distance between the two gears (m) 

Y 

displacement (m) 

Ldt 

length from the linkage 1 4 to the top surface of 




displacer (m) 

Subscripts 

Lpt 

length from the linkage 1 1 to the top surface of piston 

c 

compression chamber 


(m) 

d 

displacer 

m 

mass of working fluid (kg) 

e 

expansion chamber 

m 

mass flow rate (kg/s) 

p 

piston 

p 

pressure (kPa) 

r 

regenerative channel 

Q 

volumetric flow rate (m 3 /s) 



Q-in 

input heat transfer rate (kj/s) 

Greek symbols 

Qin 

input heat transfer per cycle (kj/cycle) 

£ 

regeneration effectiveness 

r, cp, z 

cylindrical coordinates 

y 

compression ratio 

D 

radius of displacer (m) 

V 

thermal efficiency 

r 2 

core radius of cylinder (m) 

p 

dynamic viscosity (kg/m s) 

R 

gas constant Q/kg K) 

6 

crank angle (rad) 

Rd 

offset distance from the crank to the center of gear (m) 

p 

fluid density (kg/m 3 ) 

Ra 

thermal resistance of heating head (°C/kW) 

(0 

rotational speed (rpm) 


compression space vary sinusoidally. However, it does not take into 
account non-isothermal effects and internal irreversibility caused 
by fluid friction and the pressure difference in different spaces of 
the engine. As a result, it is not capable of dealing with the ther¬ 
modynamic phenomena in a real engine. 

Recently, the present authors built a 10-W (3-type Stirling engine 
with a rhombic drive (Model 10ST1), of which the photographs are 
shown in Fig 1. According to the experiences learned from the 
development of the engine, it is recognized that major fundamental 
problems that the designers of the Stirling engines have to face 
were the time consumption in the design process and the cost of 
fabrication. In most of the cases the determination of the geomet¬ 
rical parameters of the engines is merely based on the experience of 
the designers or on insufficient costly experimental information. 
Therefore, the improvement of the engine can only be carried out in 
a relatively high-cost time-consuming trial-and-error process. In 
this regard, modeling the engines is expected to be one of the 
alternative solutions to these issues. A number of numerical models 
and simulation methods for different types of the Stirling engines 
have been proposed. The Stirling cycle machines under non- 
isothermal working condition were firstly analyzed by Finkelstein 
[9] in 1960s. His analysis was recognized as one significant devel¬ 
opment in the past several decades. His model, by means of 
introducing the heat transfer coefficient, the finite heat transfer in 
the working space was derived. The temperature variation of the 
gas in the working chambers was investigated. Schulz and 
Schwendig [10] proposed a general simulation model for Stirling 
cycles. A thermodynamic analysis of a Stirling engine including 
dead volumes of hot space, cold space and regenerator was per¬ 
formed by Kongtragool and Wongwises [11]. Most recently, Kar- 
abulut, et al. [12] presented a (3-type Stirling engine with a displacer 


driving mechanism by means of a lever. By means of nodal ther¬ 
modynamic and kinematic analysis, the instantaneous temperature 
distribution of working fluid, through the heating—cooling passage, 
conducting the cold space to hot space, was studied. Parlak, et al. 

[13] performed a thermodynamic analysis of a y-type Stirling 
engine by using a quasi steady flow model. In addition, Mahkamov 

[14] performed an axisymmetric computational fluid dynamics 
approach to the analysis of the working process of a solar Stirling 
engine, de Boer [15] tried to find the aximum attainable perfor¬ 
mance of Stirling engines and refrigerators. Andersen, et al. [16] 
presented a control-volume-based modeling in one space dimen¬ 
sion of oscillating, compressible flow in the reciprocating machines. 
Reinalter et al. [17] provided detailed performance analysis of 
a 10 kW dish/Stirling system. 

It has been recognized that these models give a certain amount 
of information that are not easily obtained by the experiments. In 
addition, modeling costs much lower since it can be readily used to 
evaluate the performance of the engine design before the design is 
submitted for fabrication. The present model is referred to as the 
second-order model. Herein the three-dimensional CFD model is 
not adopted because it is an expensive method. Though the three- 
dimensional CFD model might provide detailed information 
regarding the flow pattern and temperature and pressure distri¬ 
butions, it requires much more computation efforts than the 
present model. 

The rhombic drive, which is frequently used on a (3-type Stirling 
engine [18], utilizes ajointed rhomboid to convert linear motion of 
a reciprocating piston to a rotational of fly wheel. In the rhombic 
drive mechanism, one rigid rod is installed connecting the piston to 
the top corner of the jointed rhomboid, and another rigid rod 
connecting the displacer to the bottom corner of the rhomboid. In 
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Fig. 1. A 10-W [1-type Stirling engine with rhombic-drive (Model 10ST1) and its parts [from Power Engines and Clean Energy Laboratory (PEACE Lab.), National Cheng Rung 
University). 


addition, two symmetric gears of equal diameter are connected to 
the right and the left corners of the rhomboid fixed on the gears at 
an offset distance from gears center. The two gears are in contact 
and rotate in opposite directions. When gas pressure is applied to 
the piston, the top corner of the rhomboid is pushed downward; 
the rhomboid is flattened in the direction of the piston axis; and 
then it pushes on the gear wheels and causes them to rotate. At the 
same time, as the gear wheels rotate, the rhomboid progresses its 
change of shape and drives its bottom corner and the displacer to 
move upward. 

The purpose of the present study is aimed at development of 
a numerical model for predicting thermodynamic cycle and effi¬ 
ciency of the p-type Stirling engine with the rhombic drive, in 
which the annual gap around the displacer acts as a regenerator, so 
as to evaluate the thermodynamic performance of the engine and 
predict the periodic variation of pressures, specific volumes, 
temperatures, masses, and heat transfers in the expansion and the 
compression chambers. Note that the present model takes into 


account the non-isothermal effects, the effectiveness of the 
regenerative channel, and the thermal resistance in the heating 
head. Meanwhile, the internal irreversibility caused by pressure 
difference in the two chambers and the viscous shear effects due to 
the motion of the reciprocating displacer on the fluid flow in the 
regenerative channel gap are also included. 


2. Numerical model 

A schematic diagram of the engine is illustrated in Fig. 2. The 
working fluid contained in the engine is air. The applied model has 
been divided into three control volumes: (1) expansion chamber, 
(2) regenerative channel, and (3) compression chamber. Each of the 
three control volumes can be treated as an open system. It is 
assumed that the working fluid is an ideal gas and that the volumes 
of the chambers as well as the pressures, temperatures, and masses 
of the air in the expansion and the compression chambers are 
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varying with time transiently. However, all the time-varying 
properties are uniform inside the chambers at any time instant. 

Displacements of the piston and displacer connected with the 
rhombic-drive mechanism, Y p (t) and Y d (t), are given as: 


V e (t) = n-r 2 2 (l t -Y p (t)) (3) 

V c (t) = n-rj(Y d (t)~l d -Y p (t)) (4) 

where r 2 is the core radius of cylinder. Furthermore, differentiating 
Yp(t) and Yd(t) of Eqs. (1) and (2) with respect to t yields the 
velocities of the piston and the displacer, v p (t) and Vd(t), respec¬ 
tively, as 

i'p(t) = R d u 

x Q - 2 - KdCOS e'j 1 (5) 


"d(t) =RdO> 

x (j “ \ ~ RdC0S e ) | (6) 

2.3. Expansion chamber 

The expression of volume of the expansion chamber can be 
further written by inserting Eqs. (1)—(3) as 


cos u + sin ( 


( 3 " ( \ ~ J “ RdCOS 


-1/2 


cos t) — sin ( 


,2.A |_K dC0Sl 


-1/2 


14 = 7t rk { l t - 


L dt + ft d sin# 


-(> 3 - (k~ l 4 ~ R d cos 


1 / 2 ' 


( 7 ) 


Using the ideal-gas equation of state, the pressure in the 
expansion chamber is calculated by 


Pe 


m e RT e 


( 8 ) 


where m e is the mass and T e is the temperature of the air contained 
in the expansion chamber, and R is the gas constant.. 

Next, the air temperature in the expansion chamber (T e ) is 
calculated based on the energy conservation law: 


dU 

dt e 


Qin.e 


— 14 / 0 ut,e 



: 


O) 


Yp(t) 


L p t 3“ ivi sin 6 + 



h 

2 


R d cos 


1 1/2 


( 1 ) 


Y d (t) 


L dt + R d sind- 



U 

2 


R d cos 


1/2 


( 2 ) 


where 6 is the crank angle; i?d is the offset distance from gear 
centers to the joints between the rhomboid and the gears; / 1 , / 2 , 13 , 
and I4 represent lengths of the linkages of the rhombic-drive 
mechanism; L dt denotes the length from the linkage It to the top 
surface of the displacer; and L p[ is the length from the linkage l 4 to 
the top surface of the piston. 

The volumes of the expansion and the compression chambers, 
V e (t) and V c (t), can be calculated, respectively, in terms of Y p (t), 
Yd(f), and the cross-sectional area of the cylinder: 


where the values of h and v is chosen in accordance with the 
direction of the mass flow. Here, the mass flowing out of the 
expansion chamber is prescribed to be positive. Therefore: 

For mass leaving the expansion chamber and entering the 
regenerative channel: 

0 , h = h e , v = v e (10a) 

where h e is the specific enthalpy of the air contained in the 
expansion chamber, and v e is the average velocity of the air leaving 
the expansion chamber calculated with v e = |^p|/[p e ,r ( r 2 ~ r i)l- 
For mass entering the expansion chamber from the regenerative 
channel: 

dm 

<0. h = hj, v = Vi 

df 1 1 


(10b) 
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where hj is the specific enthalpy of the air exiting the regenerative 
channel, and Vj is the average velocity of the air calculated with 
Vi = Fafl/[Pj7t(r| - rj)]. Let 

^l e = m e c v (r e " +1 -r e ")/At (ii) 

where m e = (m" +1 + m ")/2 and At is the time step used in the 
time marching computation. In this study, At is typically assigned to 
be 1 x 1CT 6 s. 

With the help of Eqs. (10) and (11 ), Eq. (9) can be discretized into 
the finite-difference form as 


T , n+1 
1 e 


= i- 


At 

m e C v R, 


ti 

dm/. 

dF ft ' + f 


T n + At 
e m e C v 


Th 

Rn 


:(v e n+1 - v? 

At 


f dm 
for w> ° 


(12a) 


T" +1 =11- 


At 


m e C v R t 


T n + ^L 
e m e C v 


7h 

Rti 


>(V e n+1 - v e " 

At 


dm 

dt 


h. _)—L 

Hj -I- 9 


dm 

for , < 0 
dt 


(12b) 


where Th is the heat source temperature; P e = (P" +1 + P")/2; and 
the specific enthalpy h is calculated by h = C P T. Note that the input 
heat transfer rate is defined by Q ine = (r H -T e )/P t i and the 
output power due to volume expansion is determined by 
VVout.e = Pe(Vg +1 - Vq)/ At. In addition, R tl is thermal resistance of 
the heating head, and in this study R t i is set to be 0.5 °C/W. 
However, it is important to mention that the thermal resistance is 
basically dependent on the heating condition, the materials, and 
the geometry of the heating head. Hence, for a specific case, further 
experiments or numerical simulation should be required to deter¬ 
mine the value of R t i. Eqs. (12a) and (12b) are employed to calculate 
the air temperature in the expansion chamber. 


2.2. Regenerative channel 


In the present model, the gap between the outer surface of the 
displacer and the inner surface of the cylinder is referred to as the 
regenerative channel. The model of regenerator used here is 
a simplified model which assumes that the gas is incompressible 
throughout the interior chambers. The gas in the regeneration 
channel enters at one end and exits at another, and the direction 
and the mass flow rate of the gas flow in the regeneration channel 
is dependent on the reciprocating velocity of the displacer and the 
instantaneous pressure difference between the expansion and the 
compression chambers. Since the regenerative channel is usually 
rather small, it is reasonable to approximate the fluid flow in the 
channel as a fully developed flow. It implies that v r = = 0, 

i > z = vz(r). For a quasi-steady, fully developed, incompressible 
annular flow, the momentum equation can be simplified as follows. 


1 dP 
p dz 


1 9 
r dr 


_dv z 

Dr 


(13) 


Integrating the above equation yields the velocity inside the 
channel. 


vz 


r 2 dP 
4 p dz 


+ C, In r + C 2 


(14) 


where the coefficients, Ci and C 2 are determined based on the 
following boundary conditions: 

B.C-1 : r = r\, v z = i> d (t) (15a) 

B.C-2 : r = r 2 , v z = 0 (15b) 

where the reciprocating velocity of the displacer, v d (t), is deter¬ 
mined by Eq. (6). Then, the velocity profile in the regenerative 
channel is carried out as 


4,u dz 


VdW 
In (r 2 /ri) 




-In r 2 


vd(t) 

ln(r 2 /ri) 


In r 


(16) 


In the above equation, it is noticed that the effects of pressure 
difference in the two chambers and the viscous shear effects due to 
the motion of the reciprocating displacer on the fluid flow in the 
regenerative channel gap are included. The volumetric flow rate Q 
can then be calculated as 



vz(2nr) dr 


Tt / 

dP\ 

[ r 4 j-4 ( r 2 ~ r l) 2 l 

8 A 

dzy 

2 1 lnfo/ti) 


2TCr d (f) fl 
ln(r 2 /r,) [2 

v d(t) 


r 2 In r 2 - r\ In rj) - i(V 2 - r 


+ 2 7t- 


In r 2 


1 


1 


2 r2 2 ri 


ln(r 2 /n) 

The mass flow rate through the regenerative channel is 


(17) 


dm 

dt 


Pit / P e - Pc 

8 A (d 

2nv d {t)p 


ln(r 2 /ri) 


Mr2/rJ(-f n ^ lnr, )-iD 2 


+ 27t lJ d(t)P . 

ln(r 2 /r,) 


1 2 1 2 

2 r2 ~ 2 1 



(18) 


where dP/dZ has been approximated with (P e -P c )/l d . Using the 
obtained value of dm/df, the masses of the air contained in the 
expansion and the compression chambers can be updated as 


m; 


n+l 


n dm" 

= m "-dF' Af 


(19a) 


m" +1 = m" + C ^| At (19b) 

where the superscripts n and n + l represent two consecutive time 
steps. 

The mass conservation principle requires that the mass flow rate 
through the entire regenerative channel will remain constant. That 
is, 


thin.r = thout.r = th 


dm 

dt 


( 20 ) 


Applying energy conservation principle to the regenerative 
channel, one has 


Qin,r 


m(h in 


/lout) | r - th 




( 21 ) 


Due to the upwind effects, the inlet and outlet temperatures of 
the regenerative channel are different as the air flows in different 









































C.-H. Cheng, Y.-J. Yu / Renewable Energy 35 (2010) 2590-2601 


2595 


i 


t 


Expansion chamber T e 


T, 


Compression chamber T c 


Fig. 3. Inlet and outlet temperatures of air flowing through regenerative channel. (For 
the air flowing from the expansion to the compression chambers, the inlet tempera¬ 
ture is T e and outlet temperature is T,. For the air from the compression to the 
expansion chambers, the inlet temperature is T c and outlet temperature is 7}.) 



directions. Fig. 3 shows the inlet and outlet temperature of air 
flowing through the regenerative channel. For the air flowing from 
the expansion to the compression chambers, the channel inlet 
temperature is T e and channel outlet temperature is T>. For the air 
flowing from the compression to the expansion chambers, the 
channel inlet temperature is T c and channel outlet temperature is Tj. 
Note that the values of the outlet temperatures, T; and Tj, are 
strongly dependent on the effectiveness of the regenerative 
channel. 

The regenerator quality is usually defined on an enthalpy basis 
in terms of a regeneration effectiveness as follows: 


Fig. 4. Flow chart of the computation process. 

to be a function of G only. Based on the theory given by Organ [20], 
the function of e is approximated as a linear function in terms of G 
which is varied from 0.1 (at G = 0.0008 m) to 0.8 (at G = 0.0002 m). 
Hence, at G = 0.0005 m, e is calculated to be 0.45. The regeneration 
effectiveness is relatively low since the regenerator here is simply 
an annular channel. However, note that for a more complicated 
regenerator, a study of the effectiveness of the regenerator is 
definitely necessary. In consequence, one has: 

For mass flowing from the expansion to the compression 
chambers: 


actual enthalpy change through the regenerator 
maximum theoretical enthalpy change through the regenerator 


( 22 ) 


Based on the above definition of the regeneration effectiveness, 
the outlet temperatures, Tj and Tj, are determined by 

Tj = T e + e(T c - T e ), for^ > 0 (23a) 

and 

Tj = T C + e(T e - T c ), for ^ < 0 (23b) 

respectively. In general, the effectiveness of the regenerator is 
a function of the geometry and the properties of the porous 
medium used. The design of the regenerator can be optimized to 
yield a higher effectiveness [19]. As for the regenerative channel 
considered here, the regeneration effectiveness is only dependent 
on the size of the gap (G) and the channel length (1^). In this study, 
the channel length /d is fixed at 0.07946 m. Hence, e can be assumed 


> 0, T m = Te, hj n = he, Pj n = V e ; Tout = Tj, 

hout = hj, Pout = i'i (24a) 

For mass flowing from the compression to the expansion 
chambers,: 

■gg < 0, Tj n = T c , hj n = h c , Vj n = v c ; T ou t = Tj, 

h ou t = hy, Pout = Vj (24b) 

In the above equations, as already stated earlier, all the specific 
enthalpies are calculated by h = C P T based on the associated 
temperatures, and all the average velocities of the air are calculated 

by v = |^p|/[p7t(r| - r?)]. 
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Eq. (21) is used to determine the heat transfer rate (Q inr ) 
between the air flow and the walls of the regenerative channel. As 
Qtn.r > O’ ^ rneans that heat is transferred from the regenerative 
channel walls to the air, and on the contrary, as Q jnr < 0, it means 
that heat is rejected from the air to the channel walls. 

2.3. Compression chamber 

Introducing Eqs. (1) and (2) to Eq. (4) leads to the volume of the 
compression chamber: 


Table 2 

Operating variables of the base-line case (model 10ST1). 


Engine speed, to 

2000 rpm 

Compression ratio 

1.7420 

Heat source temperature, 7 h 

800 K 

Heat sink temperature, T L 

300 K 

Initial pressure, P H , Pi 

101.3 kPa 

Thermal resistance, jR fl 

0.5 °C/W 

Thermal resistance, R t 2 

0.4 °C/W 

Regenerative effectiveness, e 

0.45 


V c = Ttrl 


L dt + R d sin# - (li - ( k - l 4 - R d cos 


1/2 


Lpt + P d sin 6 + flj - R d cos 1 


2 \ 11/2 


-id 


(25) 


From the ideal-gas equation of state, the pressure in the 
compression chamber is 


Pc = 


m c RT c 

—vT 


(26) 


and, the gas temperature in the compression chamber (T c ) is 
calculated from energy equation 


d U 
dt 


a dm/, V 2 

— Qin.c ~ Wout,c + -g^r( h + ^ 


(27) 


where the values of h and v is chosen in accordance with the 
direction of the mass flow. Please be reminded that the mass 
flowing into the compression chamber is prescribed to be positive. 
Therefore: 

For mass entering the compression chamber from the regener¬ 
ative channel: 


dm 

-gp >0, h = hj, v = Vi 


(28a) 


-I 

dt < 


= m c C v (r c n+1 - T c n ) / 


At 


(29) 


where m c = (m[? +1 + m[?)/2. Note that the input heat transfer rate 
is defined by Q in c = (Tf — T c ) / R [2 . Since Ti < T c usually, the value of 
Qin, c is negative in most of the time. The output power is deter¬ 
mined by VVout.c = (P c (V c n+1 - V"))/Af. Thus, Eq. (27) reduces to 


jn+t = / 1 _ to _ 

c V rn c C v R t2 

d m(, v ? 

+ W hi + f 


T? + =^~ 

m c C v 


Ik 

R t2 


:(v ? +1 - v?) 

At 


dm 

f ° r W >0 


(30a) 


T n+1 = ( 1 - 


At 


m c C v R t 2 


T n 


At 


m c C v 


Ik 

R t2 




-Vv 


At 


dm/, v 2 

+ dT( hc + l 


f dm 
for —j— < 0 
at 


(30b) 


where Ty is the heat sink temperature; m c = (mj} +1 +m")/2 and 
P c = (P" +1 +P")/2; and R t 2 is thermal resistance of the cooling 
jacket, which is assigned to be 0.4 °C/W. Eqs. (30a) and (30b) are 
used to calculate the air temperature in the compression chamber. 

2.4. Network output and thermal efficiency 


where h, is the specific enthalpy of the air entering the compression 
chamber from the regenerative channel, and v, is the average 
velocity calculated with v t = |^f|/[Pjtt(r| - r 2 )]. 

For mass leaving the expansion chamber and entering the 
regenerative channel: 


The network output per cycle is further determined by inte¬ 
grating P with respect to V for a cycle in accordance with the P—V 
relations taken in the expansion and the compression chambers. 
That is, 


j- <0, h = h c , v = v c (28b) 

where h c is the specific enthalpy of the air contained in the 
compression chamber, and v c is the average velocity of the air 
calculated with v c = |^ ! |/[Pc ' 7T ( r 2 - r i)]- Again, the time derivative 
term on the right-hand side of equation is approximated by 


Table 1 

Geometrical variables of the base-line case (Model 10ST1). 


r 2 ( m) 

0.0205 

r, (m) 

0.02 

C(m) 

0.0005 

; I (=/ 2 = / 3 = /4)(m) 

0.018 

b (m) 

0.042 

Rd (m) 

0.0036 

L p t (m) 

0.05093 

Ldt (m) 

0.16374 

Id (m) 

0.07946 

lh (m) 

0.158 

Piston stroke (m) 

0.01 


to+t p t 0 +tp 

W out = J Pe d\4+ J P c dV c (31) 

to to 

The net heat transfer into the air per cycle is a summation of the 
heat transfers in the expansion chamber, the compression chamber, 
and the regenerative channel, which are obtained by integrating 
the heat transfer rates with respect to time for a cycle: 

Table 3 

Computation results for the base-line case (Model 10ST1). 

Work output in expansion chamber (kj/cycle) 

Work input in compression chamber (kj/cycle) 

Heat transfer in expansion chamber (kj/cycle) 

Heat transfer in compression chamber (kj/cycle) 

Heat transfer in regenerative channel (kj/cycle) 

Net heat input Qj n (kj/cycle) 

Network output W out (kj/cycle) 

Thermal efficiency 
Power output (W) 


1.236 x 1CT 3 
-7.335 x 1(T 4 
3.826 x 1(T 3 
-3.288 x 1CT 3 
-3.623 x 1CT 5 
5.025 x 10- 4 
5.025 x 10- 4 
0.131 
16.75 
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- Displac er (lower suifac e) 

-Piston (upper surface) 



time(s) 


Fig. 5. Displacement of the piston and the displacer. 


fo+fp to+tp to+tp 

Qin = J Qin,edf + J Q in ,c dt + J Q in , r dt (32) 

to to to 

where fo is an arbitrarily-selected reference time and t p represents 
the period of a cycle. Note that the energy conservation law 
requires that during the cyclic operation the network output per 
cycle should be equal to the net heat transfer into the air per cycle. 
This will be observed from the results later. 

The thermal efficiency of the Stirling engine can be expressed as 


V = 


W out 

to+tp 

J Qin.e dt 

to 


(33) 


Fig. 4 shows the flow chart of the computation process. The 
computation starts from the initial conditions for the air inside the 
engine. In the beginning, the pressure in the entire engine is set to 
be 101.3 kPa and the fluid temperature is 300 K. The computation 
then marches to the next time step and the solutions for pressures, 
masses, and temperatures in the expansion chamber, the 
compression chamber, and the regenerative channel are carried out 
by iteration at each time step. When the time exceeds the required 
maximum time, the solution process is terminated. The state of the 
system starts from the initial conditions and then achieves a stable 
periodic regime eventually. The periodic variation of the properties 
and heat transfers in different spaces of the engine can then be 
recorded and the network output and the thermal efficiency are 
further investigated. 

The geometrical and the operating variables of Model 10ST1, 
which is regarded as the base-line case in this study, are listed in 
Tables 1 and 2, respectively. In this case, the heat source and the 
heat sink temperatures are maintained at 800 and 300 K, respec¬ 
tively. The rotational speed of the engine is set at 2000 rpm, and 
the engine is operated at around atmospheric pressure. Based on 
the geometrical variables given, the compression ratio is 1.742. Note 
that some of the variables may be changed to investigate their 
effects on the performance of the engine while other variables 
are fixed. 


3. Results and discussion 

Table 3 shows the computation results for the base-line case. In 
this table, it is found that the network output per cycle 
(5.025 x 1CT 4 kj/cycle) closely agrees with the net heat transfer into 
the air per cycle. The numerical results completely comply with the 
energy conservation law. Furthermore, it is noticed that the engine 
power output is 16.75 W, and the thermal efficiency reaches merely 
13.1%. Such a poor thermal efficiency is frequently encountered 
with the small-scale engines similar to Model 10ST1. However, the 
magnitudes of the engine power output and the thermal efficiency 
may be elevated by modifying the geometrical design or by 
increasing the regeneration effectiveness of the regenerative 
channel. 

Detailed results from the present analysis of the base-line case 
are provided in Figs. 5—12. Fig. 5 illustrates the periodic displace¬ 
ment of the displacer and the piston in the cylinder for the base¬ 
line case. During the cycles, the displacer moves at a phase angle 
ahead of the piston in order to draw the working gas traversing 
back-and forth between the expansion and the compression 
chambers for heating or cooling. At some instants, the displacer and 
the piston were rather close to each other. The minimum distance 
between the piston and the displacer is treated to be an influential 
factor affecting the dead zone volume of the engine. Less minimum 
distance means smaller dead zone in the compression chamber and 
more working gas in the compression chamber being drawn into 
the expansion chamber for heating provided the piston and the 
displacer do not coincide with each other. 

Fig. 6 illustrates the volume variation of the expansion chamber 
and compression chamber. For the base-line case, the volume 
variation of the compression chamber is larger than that of the 
expansion chamber. The volume of expansion chamber is varied 
from 3 to 16.1 cm 3 , whereas the volume of compression chamber is 
from 2 to 20.3 cm 3 . In the beginning of the plot, the piston and the 
displacer are moving downwards, and the volumes of both cham¬ 
bers are increased. After a short while, the expansion chamber 
volume soon reaches its maximum and starts to decrease while the 
compression chamber volume is increased. When the volume of 
the expansion chamber reaches its minimum, the volume of the 
compression chamber reaches its maximum almost at the same 
time. 

Fig. 7 shows the pressure variation inside the expansion and 
compression chambers. The initial pressure of the two chambers is 


- Expansion chamber 

-Compression chamber 



Fig. 6. Volume variation of the expansion and the compression chambers. 
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- Expansion chamber 

-Compression chamJ 



Fig. 7. Pressure variation of the expansion and the compression chambers. 


Expansion chamber 



Fig. 9. Temperature variation in the expansion and the compression chambers. 


101.3 kPa at f=0. The pressure then is varied periodically during 
cycles. The pressure difference between these two chambers is one 
of the driving forces that make the fluid flow through the regen¬ 
erative channel. For this case, no higher pressure difference is 
observed since the engine was relatively small in size. The pressure 
difference data are shown in Fig. 8. Expansion chamber is full of 
high temperature working gas in a cycle, so the pressure inside is 
larger than that in the compression chamber usually. However, 
when the volume of the expansion chamber is further increased, 
the pressure inside is reduced to a certain low value. In that case, 
the pressure may be lower in the expansion chamber than in the 
compression chamber, and the pressure difference becomes 
negative. 

Fig. 9 illustrates the temperature variations in the expansion 
chamber and compression chamber. The heat source and the heat 
sink are maintained at 800 K and 300 K, respectively. It is inter¬ 
esting to note that the temperature in the expansion chamber could 
be higher than 800 1< at some instants, and on the contrary, the 
temperature in the compression chamber could be lower than 
300 K. This is because the air in the expansion chamber could be 
overheated due to compression and the air in the compression 


Pressure difference 



chamber could be overcooled due to expansion of the volume. This 
model applies the energy equation to the control volumes. By 
assuming the specific heat of the working gas, air, and the thermal 
resistance in the expansion chamber and compression chamber, 
the temperature variation is calculated during cycles. Since the 
temperatures of heat source and the heat sink are given, the heat 
transfer in the chambers can be calculated in terms of the thermal 
resistances and the temperature difference between the fluid 
temperature and the heat source or heat sink temperature. 

Fig. 10 shows the magnitudes of heat transfer rates in the 
expansion and the compression chambers. In this figure, the value 
of the heat transfer rate is positive indicating that the heat is 
transferred from the surroundings to the chamber. On the contrary, 
the value of the heat transfer rate is negative indicating that the 
heat is transferred from the chamber to the surroundings. It is 
observed that the heat transfer rate in the expansion chamber is not 
always positive and that in compression chamber is not always 
negative. The reason for this has already been discussed earlier. 


- Expansion chamber 

-Compression chamber 
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Fig. 11. Mass variation in the expansion and the compression chambers. 


Fig. 11 shows the variation in the masses of the working fluid in 
the expansion and compression chambers. Total amount of mass in 
the engine remains constant; therefore, as the mass in the expan¬ 
sion chamber is decreased, the mass in the compression chamber 
will be increased. 

Fig. 12 shows the P—V diagrams for the expansion chamber and 
compression chambers. The areas of the clockwise and counter¬ 
clockwise P—V cycles with the expansion and the compression 
chambers are calculated by integrating f t t ° +tp Pe dV e and 
fta +tp Pc d^c. respectively. The network output is thus obtained 
from Eq. (31 ). 

Next, effects of the size of the regenerative channel (G), heat 
source temperature (T H ), the distance between the two gears (L), 
and the offset distance from the crank to the center of gear ( Rd ) are 
studied and the results are shown in Figs. 13—16. 

Fig 13 shows the effects of the size of the regenerative gap on the 
power output and thermal efficiency of the engine. The regenera¬ 
tive gap size is varied from 0.00021 to 0.0008 m. The power output 
and thermal efficiency increase rapidly with the gap size as G is less 


-Thermal efficiency 

- Work output 



Regenerative gap (m) 

Fig. 13. Effects of the regenerative gap (G). 

than 0.0003 m. This is attributed to the significant increase in the 
shear friction loss as the gap size is reduced. On the other hand, 
when the gap size is increased, the regenerative effectiveness of the 
regenerative channels would be reduced. Therefore, there should 
be an optimal value for the gap size that leads to a peak value of the 
power output or the thermal efficiency. The curve of the power 
output indeed reaches a peak value of 16.75 Wat G = 0.0005 m. The 
thermal efficiency exhibits the same tendency, and at G = 0.0003 m 
the thermal efficiency reaches its peak value. In this figure, it is 
found that the maximum value of the thermal efficiency of the 
engine is 16.5%. 

Fig. 14 illustrates the effects of the heat source temperature (Th) 
on the engine’s power output and thermal efficiency. It is observed 
that the power output can be increased from 7.96 to 32.78 W as Th 
is elevated from 600 to 1200 K. The value of T H is strictly restricted 
by the induced thermal stress that the materials can endures. In the 
figure, the power output seems to be linearly varied with Th. 

The effects of the distance between the center of the gears (L) on 
the power output and the thermal efficiency of the engine are 
shown in Fig. 15. It is found that the power output is reduced from 



Fig. 12. P—V diagrams for the expansion and the compression chambers. 
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Fig. 14. Effects of heat source temperature (Th). 


Fig. 16. Effects of the offset distance from the crank to the center of gear (R d ). 


18.286 to 15.3 W, and the thermal efficiency from 0.1377 to 0.123, as 
the value of I is increased from 0.040 to 0.044 m. 

Fig. 16 shows the power output and the thermal efficiency as 
functions the offset distance from the crank to the center of gear 
(Rd). It is found that as Rd becomes larger, the power output and the 
thermal efficiency are both increased. Power output can be varied 
from 10.139 to 24.137 W with Rd changing from 0.0027 to 0.0045 m, 
while thermal efficiency from 0.106 to 0.150. Obviously, both the 
power output and the thermal efficiency monotonically increase 
with Rd. 

It is interesting to mention that the distance between the 
centers of the gears (L) and the offset distance from the crank to the 
center of gear (Rd) exhibit opposing effects on the performance of 
engine according to Figs. 15 and 16 These two parameters greatly 
affect the strokes of the piston and the displacer. In accordance with 
the obtained results, it is noted that the combination of these two 
parameters should be carefully designed for different requirements 
in engine performance. 


— — — - Thetmal efficiency 
- Work ot* pul 
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4. Concluding remarks 

A theoretical model which is applied to simulation of thermo¬ 
dynamic cycle and performance of a beta-type Stirling engine with 
rhombic-drive mechanism has been proposed in this study. Peri¬ 
odic variation of pressures, volumes, temperatures, masses, and 
heat transfers in the expansion and the compression chambers are 
predicted, and a parametric study of the dependence of the power 
output and thermal efficiency on the geometrical and physical 
parameters for a base-line case is performed. 

Results show that by adjusting the influential parameters 
including regenerative gap, distance between two gears, offset 
distance from the crank to the center of gear, and the heat source 
temperature, the performance of the base-line case can be 
improved. The power output of the base-line case reaches a peak 
value of 16.75 W at G = 0.0005 m, accompanied by a thermal effi¬ 
ciency of only 13.1%. If the thermal efficiency is of major concern, 
the thermal efficiency can be elevated to a peak value of 16.5% at 
G = 0.0003 m. 

It is also observed that the power output of the base-line case 
can be increased from 7.96 to 32.78 W as T H is elevated from 600 to 
1200 K. However, the value of Th is strictly restricted by the induced 
thermal stress that the materials can endures. 

Meanwhile, an increase in the distance between the center of 
the gears (L) leads to a decrease in the power output and the 
thermal efficiency. On the other hand, it is found that the power 
output and the thermal efficiency both increase with the offset 
distance from the crank to the center of gear (Rd). The two 
parameters, I and Rd, exhibit opposing effects on the performance 
of engine. 
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